Low Rank Updates for the Cholesky Decomposition

Usage of the Sherman-Morrison-Woodbury formula to update linear systems after low rank modifications of the system matrix is widespread in machine learning. However, it is well known that this formula can lead to serious instabilities in the presence of roundoff error. If the system matrix is symmetric positive definite, it is almost always possible to use a representation based on the Cholesky decomposition which renders the same results (in exact arithmetic) at the same or less operational cost, but typically is much more numerically stable. In this note, we show how the Cholesky decomposition can be updated to incorporate low rank additions or downdated for low rank subtractions. We also discuss a special case of an indefinite update of rank two. The methods discussed here are well-known in the numerical mathematics literature, and code for most of them can be found in the LINPACK suite. Note: Matlab MEX implementations for most of the techniques described here are available for download at http://www.kyb.tuebingen.mpg.de/bs/people/seeger/. If you make use of them (subject to the license), you have to cite this report and the website for obtaining the code in your publications. 1 The Problem and the Primitives Let A ∈ Rn,n be symmetric positive definite (we write A 0), i.e. xTAx > 0 for all x 6= 0. Then, A = LL for a lower triangular matrix L with positive diagonal elements, and this Cholesky decomposition is unique. Forget about ever inverting A if you don’t have to! Almost everything which you might want A−1 for can be done equally fast or faster using L, and will usually be much more numerically stable on a computer. The system Ax = b is solved as Ly = b, Lx = y which needs two back-substitutions (the procedure of solving a system with a triangular matrix).1 The operation count for a back-substitution is about half than for a matrixvector multiplication, so even if A−1 was given exactly there would be no gain in efficiency. However, the way via the Cholesky factorization in general leads to a much more accurate Confusingly (and historically) this procedure is termed differently depending on whether the system matrix is upper or lower triangular (forwardand back-substitution), although one is obtained from the other simply by inverting the vector index. We do not follow this nomenclature to avoid unnecessary confusion. solution for x. In fact, the best general way of inverting A would be to compute L and then A−1 using an algorithm based on repeated back-substitutions. Frequently occuring terms are computed as xTA−1x = ‖L−1x‖2, log |A| = 2 log |L| = 21T log(diagL). trA−1 can be computed from L in about half the time than A−1.2 Suppose now for a given statistical problem we use a representation based on a Cholesky decomposition A = LL . By a “representation” we mean that we store L explicitly in order to do back-substitutions on demand, or we maintain X = L−1B for some fixed matrix B ∈ Rn,m, or both. In the following we concentrate on this situation, although it can easily be extended to incorporate more involved usages of L. One of the simplest modifications is adding a new row/column to A. The corresponding update of L and X is fairly obvious and will not be discussed here. In this paper, we are interested in updating the representation if A is modified by adding a symmetric low rank matrix. In the simplest case, A′ = A + vvT (positive rank-1, update) or A′ = A − vvT (negative rank-1, downdate). These cases are fundamentally different in that for a positive update A′ A, as long as A is well-conditioned, our procedure will not run into trouble and A′ has larger eigenvalues than A.3 On the other hand, negative updates break down if A′ is not positive definite and can result in large errors if A′ is close to singularity. Low rank updates of the form A′ = A ± V V T , V ∈ Rn,d can always be done by applying d updates/downdates sequentially. In case of an indefinite low rank update, we do not know of a method which can be recommended in general. However, stable methods are available for special forms of indefinite rank two updates and may apply in cases of interest. Important examples are given in Section 4 and Section 5. In the remainder of this note, we assume that the factor L is given explicitly, and is to be overwritten by the new factor L′. In some applications, one is interested rather in an O(n) representation of L̃ s.t. LL̃ = L′. We comment on this issue briefly in each case, the reader will have not problem filling in the details.